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Abstract 

We investigate behavior of the continuous stochastic signals above some threshold, bursts, when the 
exponent of multiplicativity is higher than one. Earlier we have proposed a general nonlinear stochastic 
CNj model applicable for the modeling of absolute return and trading activity in financial markets which can 

be transformed into Bessel process with known first hitting (first passage) time statistics. Using these 
results we derive PDF of burst duration for the proposed model. We confirm derived analytical expressions 
by numerical evaluation and discuss bursty behavior of return in financial markets in the framework of 
modeling by nonlinear SDE. 
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1 Introduction 



Most econometric analysis of financial markets based on the various versions of stochastic differential equations 
(SDEs) have been limited to the case when the exponent of noise multiplicativity rj is lower than 1 [T]. This is 
related to the general econometric assumptions for the asset price process and existence of unique Martingale 

jjy measure 0. Nevertheless, Chan et al. [3] have shown, comparing varying econometric models of short-interest 

rate, that the models allowing multiplicativity of r\ > 1 capture the volatility changes of short-term interest 

• i-H rates better than those with n < 1. 

I We have introduced a class of non-linear stochastic differential equations (SDEs) providing time series with 

power-law statistics, and most notably reproducing 1// spectral density, 01 El E] ■ The general expression of the 
proposed class of Ito SDEs is 

dx=L- z^di. + xMW s . (1) 

qq Here x is the stochastic process exhibiting power-law statistics, r\ is the power-law exponent of the multiplicative 

noise, A is the exponent of power-law probability density function (PDF), and IF is a standard Wiener process 
(the Brownian motion) . Note that SDE was used in the modeling of trading activity and absolute return of 
—i financial markets [7] and is defined in the dimensionless scaled time, t S) we denote the scaling constant as of 

and will use the relation to the real time t of financial markets as t s = a\t. Empirical value of = 1/6 • 10 _5 (s _1 ) 
appropriate for the markets absolute return model was defined in [5J. 

Directly from the SDE ([I]) follows that the stationary probability density function (PDF) of this stochastic 
^ process is power-law, Pq(x) ~ x~ x , with the exponent A [9 . While in Refs. [TO] and later more precisely in [TT] 

it was shown that SDE provides time series with power-law spectral density 

Note that exponent of spectral density, f3, is defined only for r\ ^ 1. In case of r\ = 1 the SDE ([!]) becomes 
identical to the geometric Brownian motion and in case A = 3 one gets 1// noise. Do not be confused by 
cumulative inverse cubic law defined for return in |12j . where A = 4. 

Empirical data from the financial markets confirm the choice of equations with rj > 1 7] while modeling 
trading activity and absolute return. The case of r/ > 1 represents the higher positive signal-to-noise feedback 
and leads to the faster than exponential changes, growth or descent, observed in the time series |13) . 

Our modeling of financial market variables using the non-linear stochastic differential equations is based on 
the empirical analysis and the power-law statistics of proposed equations 7J. Providing microscopic, agent-based, 
reasoning for the proposed equations seems to be a formidable task for such complex system. Apparently, the 
development of macroscopic descriptions for the well established agent based models would be more consistent 
approach towards the understanding of the micro and macro correspondence. For such analysis one should 
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Figure 1: Time series exhibiting bursty behavior, I(t). Here hj is threshold value, above which bursts are 
detected, ti is the three visible threshold passage events, I m ax is the highlighted burst's peak value. The other 
relevant statistical properties are defined in respect to those shown in this figure: T = ti — t\, 6 = t% — t%, 
t = T + 6 =t 3 -*i- 

select simple agent based models with established or expected macroscopic description. First of all we expect 
that herding, contagion and cascading behavior of agents should lead to the macroscopic description with the 
non-linear stochastic differential equations. 

The convenience of this approach is already confirmed by some recent publications considering this problem 
|14i ri5, 6J. Kirman's ant colony model |16j is an agent-based model, which explains the importance of herding. 
As human crowd behavior ideologically is very similar, the ant colony model actually was built as a general 
framework for financial market modeling jTTJ [T51 [TB]. In [TH IT51 [5] we follow the works by Alfarano et al. 
[TTl ITS] to introduce model modifications, which allow us to obtain agent-based models for the financial markets 
exhibiting a macroscopic description in terms of SDEs with rj > 1. 

In this contribution we demonstrate that the general class of SDEs ([!]) can be transformed to the Bessel 
process, which represents a special family of diffusion models applicable in econometric analysis pQ. This allows 
us to derive explicit form of burst statistics generated by the SDE ([!]) . Exponent of multiplicativity rj is a key 
parameter of defined statistics. We also present analysis of empirical data providing an evidence that burst 
statistics of return in financial markets can be modeled by non-linear stochastic differential equations with 
multiplicativity as high as r/ = | . 

2 Statistics of stochastic bursty time series 

Looking for the best form of SDE applicable to the modeling of complex systems one has to use the most distinct 
features of the observed behavior. Spikes, bubbles or bursts of the observed signals is one of the characteristic 
features of complex systems, as well as 1/f noise. Well defined statistics of bursty behavior can serve as an 
additional source of information about the system. It is natural to expect that parameter of multiplicativity 
rj in Eq. (JlJ is responsible for the bursty behavior of defined stochastic signal. Therefore described statistics 
of bursty behavior would help establish relation between model parameters and empirical data. Although we 
concentrate here more on the return not price dynamics, such characteristics of bursts as duration and peak 
value should be indispensable in risk valuation process. 

We define a burst as a part of time series lying above certain threshold, h. In Fig. [T]we have presented an 
example burst of the simple bursty time series, I(t) (hence the threshold is denoted as hi). Evidently a burst 
itself can be characterized by its duration, T = t 2 — t±, peak value, I m axi and burst size, defined as the area 
above the threshold yet bellow time series (highlighted), S. One can also introduce inter-burst, 9 = £3 — 12, and 
waiting, r = t% — t\, times to be able to fully grasp the statistical features of bursty behavior. We will derive 
the explicit form PDF of burst duration T, generated by Eq. ([!]) , and will analyze numerically correspondence 
with empirical data of returns and relations to the other statistical characteristics of bursty behavior. 

There is a well established hitting (passage) time framework, which is of high interest in both mathematical 
finance [T] and physics (SJ[TS]. This framework can be also applied in the description of burst duration PDF. 
Actually first hitting time of the stochastic process starting infinitesimally near above the hitting threshold is 
the same as burst duration. It is so as the first passage of threshold ends the burst by its definition. 
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2.1 Numerical and empirical definitions of bursty behaviour statistics 

In this contribution we consider stochastic model driven by SDE (JlJ and empirical time series of absolute return. 
Let us present a brief discussion on how we deal with the numerical and empirical time series. 

SDE ([!]) serves as a most simple definition of very wide class of stochastic processes with power law statistics 
for high values of the signal intensity, x. Any real system with corresponding power law statistics of high values 
has to be restricted from the side of small values. One can consider reflecting boundary condition at the point 
x = as the most simple case. Nevertheless, the more general approach can be implemented by wide choice of 
compressed exponential restrictions, which can be introduced by additional term in SDE 

dx = (v - ± + fff) x 2 ^dt s + zW s . (3) 

Here x m i n denotes the diffusion restriction point and m defines how sharp restriction is. Steady state distribution 
of x is as follows 



mxt 1 



n i 



As we analyze bursts, behavior in the region of high values of x, it is reasonable to assume x m i n = 1 and 
m = 2. With such assumption we numerically solve ^ by evaluating the difference equations: 

x l+ i = Xi + k 2 \ r\ - \ + j^j Xi + Ky/xid, (5) 

tsA+l = U + x 2r,-2 , (6) 

here Q is a normally distributed random variable with unit variance and zero mean, n is a numerical precision 
parameter (it should be significantly smaller than 1). The above difference equations follow from the Euler- 
Maruyama method for the numerical solution of the stochastic differential equations |20j with variable time 
steps. Looking for the higher precision one may alternatively choose Milstein method |20j . 

We iterate through the ^ and Q until we obtain time series with set amount of bursts. The suitable 
precision was achieved from time series with 10 5 bursts. We choose threshold value h x — 2 to ensure sufficient 
number of burst and A = 4 as well defined empirical exponent of return distribution |12j , in our all numerical 
calculations. 

Our empirical data set includes all trades made on NYSE, which were made from January, 2005 to March, 
2007 and involved 24 different stocks, ABT, ADM, BMY, C, CVX, DOW, FNM, GE, GM, HD, IBM, JNJ, 
JPM, KO, LLY, MMM, MO, MOT, MRK, SLE, PFE, T, WMT, XOM. We have shown in [7 that the more 
sophisticated versions of (§ may be well used to model absolute return and trading activity of different stocks 
from NYSE and Vilnius Stock Exchange. In this approach the normalized absolute return time series of different 
stocks served as independent realizations of the stochastic process driven by the same nonlinear SDE with the 
same set of parameters. Such universal nature of stock market return behavior is hidden under secondary high 
frequency stochastic process with the g-Gaussian PDF. 

The double stochastic nature of the return makes analysis of the bursty behavior in empirical stock data 
as a more complex task. One minute return time series are very noisy, see [7] for more details. This noise 
of normalized time series can be diminished by using moving average filter. We have used one hour window 
moving average filter on the empirical data of one minute absolute returns. Time series of all stocks were used 
to define statistics of burst durations, peak values and sizes. 



2.2 Obtaining probability density function of burst durations using the first hit- 
ting times of Bessel process 

There are few simple and highly applicable models for which hitting times statistics are known. These models 
include, but are not limited to, Brownian motion, geometric Brownian motion and Bessel process [T]. The 
Bessel process, 

is one of the most interesting as some prominent mathematical finance models can be transformed to a similar 
form. Bessel process is also interesting in Physics since it describes the evolution of the Euclidean norm, R — \R\, 
of the Y-dimensional Brownian particle. The widely used parameter v, 

N 

" = -2-1. (8) 

is known as the index of the Bessel process related to the TV-dimensional Brownian diffusion. Note that for 
N > 1, or alternatively v > —0.5, R tends to diverge towards infinity. 
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One can reduce ([I]) to Bessel process by using Lamperti transformation, i : x i— > y(x). The exact form of 
transformation can be obtained by requiring that 



x^d x y(x) = ±1. (9) 

Note that this requirement follows from the the Ito variable substitution formula 9j. Thus the Lamperti 
transformation has the following form 

I:i4i/(i) = 7 ^ T . (10) 

yy ' (rj- l)^- 1 v ; 

In such case the stochastic differential equation ([!]) can be transformed to the Bessel process, 

dy=(u+l^ + dW s , (11) 

with index v = . The dimension of the obtained Bessel process is given by N = 2(v + 1) = ^Ej- 

Let us assume that a burst starts at time to, thus xq = x(to) exceeds the threshold, h x , by some small 
amount. The burst lasts until x(t) crosses back h x , from the above. Equivalently, in the terms of Bessel process 
the burst lasts until at a certain time, t, the y process crosses the boundary h y := £(h x ) from the below, while 
the starting position, yo = y(to), hi terms of Bessel process also lies below the threshold, yo — £(xq) < h y . 
We see that the burst durations of the stochastic process driven by SDE ([!]) is related to the well known first 
passage times of the Bessel process. 

By choosing y arbitrarily close yet below h y , we can obtain an estimate for the burst duration, T, in terms 

of first passage times of the Bessel process, : 



T = r M 

yo,h y 



inf it, y(t) >h y \, 0<h y -y Q <e, (12) 
t>to v. ) 

where e is an arbitrary small positive constant. As given in |21j . the following holds for < yo < h y : 

Ay) (4-\ _ \ \ n y J 

y ° k=l 

where p y ^ h (t) is a probability density function of the first passage times at level h y of Bessel process with 
index v starting from yo, J v is a Bessel function of the first kind of order v and j Uj k is a k-th zero of J„. 

We have to replace p^ a h (t) by density function regarding h y to avoid trivial convergence of p y l> J fl {t) to zero, 
when yo — > h y . This is achieved introducing PDF (t) as a probability density function of burst duration 




Phy W= T — » 

* va^h y tl y — yo 

where the threshold is set at level h y and v is related to the model parameters (see discussion on the transfor- 
mation to Bessel process). To evaluate the limit we have to expand J„ \J^-j v ,kJ near J^- = 1: 

Jv (j^jv,k \ ~ Jv{jv,k) ~ [}~ h~) \ v Jv{jv,k) - jv,kJl+v{jv,k)] = 

- (1 - £)3u,kJi+„Uu,k)- (15) 



By using this expansion we can rewrite (141 as: 

oo 

< } (*) = criX;^«p(-Sf*)' (16) 



k=l 




here C\ is a normalization constant. Since are almost equally spaced [22], we can replace the sum by 
integration 
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Figure 2: Simple SDE (a) and empirical versus the double stochastic model (b) PDF of the burst durations, 
h x = 2. In sub-figure (a) numerical data is represented by the filled shapes, while fits, (17), are represented by 



gray curves. Simple SDE, parameters were set as follows: erf = 1/6 • 10 s (in all three cases), A = 4 (in 
all three cases), 77 = 2.5 (squares, v = 0), rj = 2 (circles, v — 0.5) and r\ — 1.5 (triangles, v = 2). Empirical data 
(gray empty circles in the sub- figure (b) is plotted versus the numerical data from the complex SDE, ( 20 ), (black 
empty squares) and from the double stochastic model (black filled squares). Complex SDE, (20), parameters 
were set as follows: of = 1/6 ■ 10 _5 (s _1 ), rj = 2.5, A = 3.6, e = 0.017, x max = 10 3 . Double stochastic model 
uses the same parameters as complex SDE and the additional parameters, which were set as follows: fo = 0.4, 
A 2 = 5. 



Note that derived PDF in the form of Eqs. ( 16 1 or ( 17 ) diverge, when t approaches zero. Therefore normalization 
constants C\ and C2 can be defined if some minimum value of t is supposed. Accuracy of numerical calculations 
or minimum intertrade time can be considered as possible choices. From the above follows that the distribution 
of burst durations of SDE ([!]) can be approximated by a power-law with exponential cut-off according to 

pS } (*)~*~ 3/2 . whcn (18) 



<V) when *»g (is) 

This result is in agreement with a general property of one dimensional diffusion processes presented in |19j . 
namely that the asymptotic behavior of first hitting times is a power-law i~ 3 / 2 irrespectively of the nature 
of stochastic process or the actual form of Langevin and Fokker-Plank equations. The exponential cutoff for 
longer burst durations is caused by the direction preference, note the positive drift term in case of N > 1, or 
alternatively v > —0.5, of Bessel processes. Numerical solutions of the SDE ([!]) confirm the derived probability 
density function, (17), of the burst duration, T, see Fig. [2] (a). 

Empirical data, as shown in Fig. [2] (b), also has similar asymptotic behavior for short and long burst 
durations, though the fitting using ( |17[ ) would be inconsistent for the intermediate burst durations, note the 
cusp. There are a numerous reasons for this. Firstly, we were unable to remove intra-day pattern from the 
empirical time series. But the main reason is that in order to reproduce the correct shape of the empirical 
probability density function for the intermediate values of T one must use the double stochastic model, driven 
by a more sophisticated version of the SDE ([!]) . In Fig. [2] (b) we demonstrate pretty good agreement of the 
probability density function of T retrieved from the empirical NYSE time series and a more sophisticated double 
stochastic return model discussed in Appendix A. 

2.3 Power law interdependencies between burst related variables 

Another interesting feature of both numerical results and empirical data is that burst duration, peak value and 
burst size are correlated or even interdependent, (T, x max , S). The scatter plots of these variables also reveal 
power-law asymptotic behavior - x max cx T$ (see Fig. |3J) and S oc Ts (see Fig. Ul). Consequently, as follows 



cxp 



2h: 



from these relations S cx (xmax)% oc Xmax (see Fig. [5]). From SDE |l| numerically defined asymptotic relations 
are shown to hold for both empirical data and numerical results, though the theoretical power-law holds only 
in narrow region of empirical data. The interdependence between these variables suggest that the geometry 
of burst remains qualitatively the same in the fitted (T, x max ,S) regions. Theoretical consideration of these 
relations in the region of T PDF exponential cut-of will be continued. 
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Figure 3: The scatter plot of the numerical (a) and empirical (b) burst peak values vs burst durations, h x = 2. 
In both sub-figures filled squares represent mean values and error bars - variability, while power-law fits, a = |, 
are represented by gray curves. Model, ([l]), parameters were set as follows: of = 1/6 • 10 _5 s _1 , rj — 2. A = 4 
(i/ = Q). 
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Figure 4: The scatter plot of the numerical (a) and empirical (b) burst size vs burst durations, h x = 2. In both 
sub-figures filled squares represent mean values and error bars - variability, while power-law fits, a = |, are 
represented by gray curves. Model, ([I]), parameters are the same as in Fig. [3] 




Figure 5: The scatter plot of the numerical (a) and empirical (b) burst size vs burst peak values, h x = 2. In 
both sub-figures filled squares represent mean values and error bars - variability, while power-law fits, a = § 
are represented by gray curves. Model, 0, parameters are the same as in Fig. [3] 
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3 Conclusions 



Starting from the general expression of nonlinear SDE ([I]) and well-known PDF expression of first hitting time 



for Bessel process we have derived PDF of burst duration (16) and its approximation ([17]) for the nonlinear 
stochastic process driven by ([I]). The understanding of burst statistics is needed to enable efficient modeling 
of trading activity and return in financial markets based on the various versions of SDE [7 . We do expect 
that proposed class of SDE is applicable in the modeling of other complex systems and defined statistics of 
bursts can be helpful in risk analyzes. We have also confirmed analytical results of burst statistics by numerical 
calculations of SDE ([!]), and compared it to the empirical data of return in financial markets. The obtained 
results encourage us to continue research of burst statistics in the financial markets seeking to develop more 
sophisticated versions of stochastic models. Most importantly we find that burst related statistical properties 
are related to the key parameter - exponent of stochastic multiplicativity r\ characterizing the dynamics of risk 
in financial markets. Both numerical results and empirical data also reveal power-law asymptotic behavior of 
burst duration, peak value and burst size scatter plots. The exponents of observed power-law statistics probably 
have universal nature and will be considered in our future work. 

Our analysis of empirical data suggests that long range fluctuations of absolute return in financial markets 
can be modeled by non-linear stochastic differential equations with 77 up to |. 
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A Double stochastic return model 



The class of equations based on SDE |T]) gives only a general idea how to model power-law statistics of trading 
activity and return in the financial markets. The problem is to determine the parameter set A and 77, which 
would enable reproduction of the empirical values of A and j3. The task becomes even more complicated if 
one considers the more sophisticated trends in the spectral density - power spectral densities have not a single, 
but two power-law regions with different values of /3. In the series of papers [331 H] we have shown that 
trading activity and return can be modeled by a more sophisticated stochastic differential equation than ([!]) 
now including the two powers of the noise multiplicativity. In the case of return instead of Eq. ([I]) one should 
use [5] 



da; 



77 



A 
2 



(1 + x 2 y> 



-xd/j. 



(1 + x 2 



-dW s 



(20) 



(e-v/l + x 2 + l) 2 eVl + x 2 + l 

here e divides the area of diffusion into the two different noise multiplicativity regions to ensure the spectral 

density of |a;| with two power law exponents, term (^ x x ^ gives the exponential diffusion restriction for large 

values of variable when x > x max . 

The proposed form of the more complex SDE enables reproduction of the main statistical properties of the 
return observed in the financial markets. This provides an approach to the financial markets with behavior 
dependent on the level of activity and exhibiting two stages: calm and excited. One more peculiarity of 
the proposed model is that signal x serves only as modulating one of secondary high frequency (/-Gaussian 
fluctuations. We formalized the empirical return r t in the model as instantaneous q-Gaussian fluctuations £ 
with a slowly diffusing parameter r and constant A = 5 



n = £{n>, A}. 



q-Gaussian distribution of £ can be expressed as follows: 



P ro ,x(r) = 



(21) 



(22) 



The variable ro serves as a measure of instantaneous volatility of high frequency return fluctuations |S]. We 
proposed to model the measure of volatility ro by the scaled continuous stochastic variable x, which can be 
interpreted as average return per unit time interval. Through the empirical analysis of high frequency empirical 
data from NYSE [8] we have introduced the following relation 



r (t,r) = 1 + 



t B +T B 



x(s)ds 



(23) 



7 



where tq is an empirical parameter and the average return per unit time interval x(t s ) can be modeled by a 
nonlinear SDE (20), expressed in scaled dimensionless time t s = aft. 
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